library(rethinking)
Loading required package: rstan
Loading required package: StanHeaders
Loading required package: ggplot2
rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Loading required package: parallel
rethinking (Version 2.11)

Attaching package: ‘rethinking’

The following object is masked from ‘package:stats’:

    rstudent
library(brms)
Loading required package: Rcpp
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Loading 'brms' package (version 2.13.5). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: ‘brms’

The following objects are masked from ‘package:rethinking’:

    LOO, stancode, WAIC

The following object is masked from ‘package:rstan’:

    loo

The following object is masked from ‘package:stats’:

    ar
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ tibble  3.0.3     ✓ dplyr   1.0.0
✓ tidyr   1.1.0     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0
✓ purrr   0.3.4     
── Conflicts ──────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x tidyr::extract() masks rstan::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x purrr::map()     masks rethinking::map()
germ <- read_csv("../Dimensions/hydrothermal-all-species/data/light_round1_tall.csv") %>% filter(wps==0) %>%
Warning message:
Unknown or uninitialised column: `germ`. 
  select(pops, temps, total_seeds, germ, day, cumulative_germ)
Parsed with column specification:
cols(
  pops = col_character(),
  temps = col_double(),
  wps = col_double(),
  date = col_character(),
  total_seeds = col_double(),
  germ = col_double(),
  start_date = col_date(format = ""),
  census_date = col_date(format = ""),
  day = col_double(),
  cumulative_germ = col_double(),
  cumulative_prop_germ = col_double()
)
germ

what if need one event per row:

one_per_row <- function(df) {
  total_seed <- max(df$total_seeds, sum(df$germ))
  newdata <- tibble(id=1:total_seed, germ=0, day=max(df$day))
  df <- df %>% filter(germ>0)
  count <- 1
  if (nrow(df) > 0) {
    for (i in 1:nrow(df)) { # we look at each row of the df where germination occured
      for (j in 1:df$germ[i]) { # now update the newdata to reflect the germiantion of each seed
        newdata$germ[count] <- 1
        newdata$day[count]=df$day[i]
        count <- count+1 # count keeps track of which individual we are at in the new data
      } # for j
    } # for i
  } # if 
  return(newdata)
}

germone <- germ %>% group_by(pops, temps) %>%
  select(-cumulative_germ) %>% # not needed in this encoding (I think...in any case would need to be recalculated)
  nest() %>%
  mutate(newdata=map(data, one_per_row)) %>%
  select(-data) %>%
  unnest(newdata)

germone

##1

Is effect of temp ~linear on cum germ?

germ %>% filter(pops=="STDI", day==28) %>%
Warning message:
Unknown or uninitialised column: `germ`. 
  ggplot(aes(x=temps,y=cumulative_germ)) +
  geom_col()

no, so we can’t treat temp as a continuous linear predictor

germ.stdi <- germone %>% filter(pops=="STDI") %>% select(-pops)
Adding missing grouping variables: `pops`
germ.stdi

Try a censored lambda model

d <- list(germ=germ.stdi$germ, temps=as.numeric(as.factor(germ.stdi$temps)), day=germ.stdi$day)
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
3: Unknown or uninitialised column: `germ`. 
m1.1 <- ulam(
  alist(
    day | germ==1 ~ exponential( lambda),
    day | germ==0 ~ custom(exponential_lccdf( !Y | lambda)),
    lambda <- 1.0 / mu,
    log(mu) <- a[temps],
    a[temps] ~ normal(0,1)),
  data=d,
  chains=4,
  cores = 4
)
starting worker pid=54244 on localhost:11220 at 11:44:48.270
starting worker pid=54258 on localhost:11220 at 11:44:48.609
starting worker pid=54273 on localhost:11220 at 11:44:49.300
starting worker pid=54288 on localhost:11220 at 11:44:49.607

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000109 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.09 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 7.5e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 7.4e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.74 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.294589 seconds (Warm-up)
Chain 1:                0.230801 seconds (Sampling)
Chain 1:                0.52539 seconds (Total)
Chain 1: 
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)

SAMPLING FOR MODEL '6c7ed290b113c205c3adee6a5a6dddff' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 6.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.280862 seconds (Warm-up)
Chain 2:                0.212282 seconds (Sampling)
Chain 2:                0.493144 seconds (Total)
Chain 2: 
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.310661 seconds (Warm-up)
Chain 3:                0.219301 seconds (Sampling)
Chain 3:                0.529962 seconds (Total)
Chain 3: 
Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.323819 seconds (Warm-up)
Chain 4:                0.215416 seconds (Sampling)
Chain 4:                0.539235 seconds (Total)
Chain 4: 
precis(m1.1, depth = 2)
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 

posterior

preddata <- expand_grid(temps=1:8, realtemps=sort(unique(germ.stdi$temps)), day=1:28)
Warning message:
Unknown or uninitialised column: `germ`. 
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
 $ lambda: num [1:2000, 1:1792] 0.00439 0.00507 0.00342 0.00459 0.00509 ...
 $ mu    : num [1:2000, 1:1792] 228 197 293 218 197 ...

mu is average days to germ, lambda is rate parameter. neither change over time, of course, so having days doesn’t really make sense.

preddata$mu <- apply(pred$mu,2,mean)
Warning message:
Unknown or uninitialised column: `germ`. 
preddata$low <- apply(pred$mu,2,HPDI)[1,]
preddata$high <- apply(pred$mu, 2, HPDI)[2,]

Single temp. Don’t need day posterior

preddata <- expand_grid(temps=3)
Warning message:
Unknown or uninitialised column: `germ`. 
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
 $ lambda: num [1:2000, 1] 0.0176 0.0212 0.0178 0.0241 0.0221 ...
 $ mu    : num [1:2000, 1] 56.9 47.3 56 41.5 45.3 ...

how to convert to probs? use pexp.

predprobs <- pexp(1:28,rate=pred$lambda[1])
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
plot(x=1:28,y=predprobs, type="l")

even thought it isn’t using day, including day in the prediction data frame will help me keep data in the correct format. Maybe.

preddata <- expand_grid(temps=1:8, day=1:28)
Warning message:
Unknown or uninitialised column: `germ`. 
pred <- link(m1.1, data = preddata)
str(pred)
List of 2
 $ lambda: num [1:2000, 1:224] 0.00439 0.00507 0.00342 0.00459 0.00509 ...
 $ mu    : num [1:2000, 1:224] 228 197 293 218 197 ...
predresults <- preddata %>%
  mutate(lambda=as.list(as.data.frame(pred$lambda)))
predresults
predresults <- predresults %>%
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
  mutate(probs=map2(day, lambda, ~ pexp(.x, .y)),
         mu=map_dbl(probs, mean),
         lower=map_dbl(probs, ~ HPDI(.)[1] %>% unlist()),
         upper=map_dbl(probs, ~ HPDI(.)[2]) %>% unlist())
predresults
predresults %>% select(-lambda, -probs) %>%
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
3: Unknown or uninitialised column: `germ`. 
4: Unknown or uninitialised column: `germ`. 
5: Unknown or uninitialised column: `germ`. 
6: Unknown or uninitialised column: `germ`. 
7: Unknown or uninitialised column: `germ`. 
8: Unknown or uninitialised column: `germ`. 
9: Unknown or uninitialised column: `germ`. 
  mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
  ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
  geom_line() 

Add realdata:

stdi.plot <- germ %>% filter(pops=="STDI") %>% 
Warning messages:
1: Unknown or uninitialised column: `germ`. 
2: Unknown or uninitialised column: `germ`. 
3: Unknown or uninitialised column: `germ`. 
  select(day, temps, cumulative_germ, total_seeds) %>%
  mutate(temps=as.factor(temps),
         prop.germ=cumulative_germ/total_seeds)

predresults %>% select(-lambda, -probs) %>%
  mutate(temps=factor(temps, labels=as.character(sort(unique(germ.stdi$temps))))) %>%
  ggplot(aes(x=day,y=mu,color=temps,group=temps)) +
  geom_line() +
  geom_point(aes(y=prop.germ), data=stdi.plot)

Poor!

LS0tCnRpdGxlOiAiZ2VybWluYXRpb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyfQpsaWJyYXJ5KHJldGhpbmtpbmcpCmxpYnJhcnkoYnJtcykKbGlicmFyeSh0aWR5dmVyc2UpCmBgYAoKYGBge3J9Cmdlcm0gPC0gcmVhZF9jc3YoIi4uL0RpbWVuc2lvbnMvaHlkcm90aGVybWFsLWFsbC1zcGVjaWVzL2RhdGEvbGlnaHRfcm91bmQxX3RhbGwuY3N2IikgJT4lIGZpbHRlcih3cHM9PTApICU+JQogIHNlbGVjdChwb3BzLCB0ZW1wcywgdG90YWxfc2VlZHMsIGdlcm0sIGRheSwgY3VtdWxhdGl2ZV9nZXJtKQpnZXJtCmBgYAoKd2hhdCBpZiBuZWVkIG9uZSBldmVudCBwZXIgcm93OgoKYGBge3J9Cm9uZV9wZXJfcm93IDwtIGZ1bmN0aW9uKGRmKSB7CiAgdG90YWxfc2VlZCA8LSBtYXgoZGYkdG90YWxfc2VlZHMsIHN1bShkZiRnZXJtKSkKICBuZXdkYXRhIDwtIHRpYmJsZShpZD0xOnRvdGFsX3NlZWQsIGdlcm09MCwgZGF5PW1heChkZiRkYXkpKQogIGRmIDwtIGRmICU+JSBmaWx0ZXIoZ2VybT4wKQogIGNvdW50IDwtIDEKICBpZiAobnJvdyhkZikgPiAwKSB7CiAgICBmb3IgKGkgaW4gMTpucm93KGRmKSkgeyAjIHdlIGxvb2sgYXQgZWFjaCByb3cgb2YgdGhlIGRmIHdoZXJlIGdlcm1pbmF0aW9uIG9jY3VyZWQKICAgICAgZm9yIChqIGluIDE6ZGYkZ2VybVtpXSkgeyAjIG5vdyB1cGRhdGUgdGhlIG5ld2RhdGEgdG8gcmVmbGVjdCB0aGUgZ2VybWlhbnRpb24gb2YgZWFjaCBzZWVkCiAgICAgICAgbmV3ZGF0YSRnZXJtW2NvdW50XSA8LSAxCiAgICAgICAgbmV3ZGF0YSRkYXlbY291bnRdPWRmJGRheVtpXQogICAgICAgIGNvdW50IDwtIGNvdW50KzEgIyBjb3VudCBrZWVwcyB0cmFjayBvZiB3aGljaCBpbmRpdmlkdWFsIHdlIGFyZSBhdCBpbiB0aGUgbmV3IGRhdGEKICAgICAgfSAjIGZvciBqCiAgICB9ICMgZm9yIGkKICB9ICMgaWYgCiAgcmV0dXJuKG5ld2RhdGEpCn0KCmdlcm1vbmUgPC0gZ2VybSAlPiUgZ3JvdXBfYnkocG9wcywgdGVtcHMpICU+JQogIHNlbGVjdCgtY3VtdWxhdGl2ZV9nZXJtKSAlPiUgIyBub3QgbmVlZGVkIGluIHRoaXMgZW5jb2RpbmcgKEkgdGhpbmsuLi5pbiBhbnkgY2FzZSB3b3VsZCBuZWVkIHRvIGJlIHJlY2FsY3VsYXRlZCkKICBuZXN0KCkgJT4lCiAgbXV0YXRlKG5ld2RhdGE9bWFwKGRhdGEsIG9uZV9wZXJfcm93KSkgJT4lCiAgc2VsZWN0KC1kYXRhKSAlPiUKICB1bm5lc3QobmV3ZGF0YSkKCmdlcm1vbmUKYGBgCgojIzEKCklzIGVmZmVjdCBvZiB0ZW1wIH5saW5lYXIgb24gY3VtIGdlcm0/CgpgYGB7cn0KZ2VybSAlPiUgZmlsdGVyKHBvcHM9PSJTVERJIiwgZGF5PT0yOCkgJT4lCiAgZ2dwbG90KGFlcyh4PXRlbXBzLHk9Y3VtdWxhdGl2ZV9nZXJtKSkgKwogIGdlb21fY29sKCkKYGBgCm5vLCBzbyB3ZSBjYW4ndCB0cmVhdCB0ZW1wIGFzIGEgY29udGludW91cyBsaW5lYXIgcHJlZGljdG9yCgpgYGB7cn0KZ2VybS5zdGRpIDwtIGdlcm1vbmUgJT4lIGZpbHRlcihwb3BzPT0iU1RESSIpICU+JSBzZWxlY3QoLXBvcHMpCmdlcm0uc3RkaQpgYGAKVHJ5IGEgY2Vuc29yZWQgbGFtYmRhIG1vZGVsCmBgYHtyfQpkIDwtIGxpc3QoZ2VybT1nZXJtLnN0ZGkkZ2VybSwgdGVtcHM9YXMubnVtZXJpYyhhcy5mYWN0b3IoZ2VybS5zdGRpJHRlbXBzKSksIGRheT1nZXJtLnN0ZGkkZGF5KQptMS4xIDwtIHVsYW0oCiAgYWxpc3QoCiAgICBkYXkgfCBnZXJtPT0xIH4gZXhwb25lbnRpYWwoIGxhbWJkYSksCiAgICBkYXkgfCBnZXJtPT0wIH4gY3VzdG9tKGV4cG9uZW50aWFsX2xjY2RmKCAhWSB8IGxhbWJkYSkpLAogICAgbGFtYmRhIDwtIDEuMCAvIG11LAogICAgbG9nKG11KSA8LSBhW3RlbXBzXSwKICAgIGFbdGVtcHNdIH4gbm9ybWFsKDAsMSkpLAogIGRhdGE9ZCwKICBjaGFpbnM9NCwKICBjb3JlcyA9IDQKKQpgYGAKCmBgYHtyfQpwcmVjaXMobTEuMSwgZGVwdGggPSAyKQpgYGAKCnBvc3RlcmlvcgpgYGB7cn0KcHJlZGRhdGEgPC0gZXhwYW5kX2dyaWQodGVtcHM9MTo4LCBkYXk9MToyOCkKcHJlZCA8LSBsaW5rKG0xLjEsIGRhdGEgPSBwcmVkZGF0YSkKc3RyKHByZWQpCmBgYAoKbXUgaXMgYXZlcmFnZSBkYXlzIHRvIGdlcm0sIGxhbWJkYSBpcyByYXRlIHBhcmFtZXRlci4gIG5laXRoZXIgY2hhbmdlIG92ZXIgdGltZSwgb2YgY291cnNlLCBzbyBoYXZpbmcgZGF5cyBkb2Vzbid0IHJlYWxseSBtYWtlIHNlbnNlLgoKYGBge3J9CnByZWRkYXRhJG11IDwtIGFwcGx5KHByZWQkbXUsMixtZWFuKQpwcmVkZGF0YSRsb3cgPC0gYXBwbHkocHJlZCRtdSwyLEhQREkpWzEsXQpwcmVkZGF0YSRoaWdoIDwtIGFwcGx5KHByZWQkbXUsIDIsIEhQREkpWzIsXQpgYGAKClNpbmdsZSB0ZW1wLiAgRG9uJ3QgbmVlZCBkYXkKcG9zdGVyaW9yCmBgYHtyfQpwcmVkZGF0YSA8LSBleHBhbmRfZ3JpZCh0ZW1wcz0zKQpwcmVkIDwtIGxpbmsobTEuMSwgZGF0YSA9IHByZWRkYXRhKQpzdHIocHJlZCkKYGBgCgpob3cgdG8gY29udmVydCB0byBwcm9icz8gdXNlIHBleHAuCgpgYGB7cn0KcHJlZHByb2JzIDwtIHBleHAoMToyOCxyYXRlPXByZWQkbGFtYmRhWzFdKQpwbG90KHg9MToyOCx5PXByZWRwcm9icywgdHlwZT0ibCIpCmBgYAoKZXZlbiB0aG91Z2h0IGl0IGlzbid0IHVzaW5nIGRheSwgaW5jbHVkaW5nIGRheSBpbiB0aGUgcHJlZGljdGlvbiBkYXRhIGZyYW1lIHdpbGwgaGVscCBtZSBrZWVwIGRhdGEgaW4gdGhlIGNvcnJlY3QgZm9ybWF0LiAgTWF5YmUuCgpgYGB7cn0KcHJlZGRhdGEgPC0gZXhwYW5kX2dyaWQodGVtcHM9MTo4LCBkYXk9MToyOCkKcHJlZCA8LSBsaW5rKG0xLjEsIGRhdGEgPSBwcmVkZGF0YSkKc3RyKHByZWQpCmBgYAoKYGBge3J9CnByZWRyZXN1bHRzIDwtIHByZWRkYXRhICU+JQogIG11dGF0ZShsYW1iZGE9YXMubGlzdChhcy5kYXRhLmZyYW1lKHByZWQkbGFtYmRhKSkpCnByZWRyZXN1bHRzCmBgYApgYGB7cn0KcHJlZHJlc3VsdHMgPC0gcHJlZHJlc3VsdHMgJT4lCiAgbXV0YXRlKHByb2JzPW1hcDIoZGF5LCBsYW1iZGEsIH4gcGV4cCgueCwgLnkpKSwKICAgICAgICAgbXU9bWFwX2RibChwcm9icywgbWVhbiksCiAgICAgICAgIGxvd2VyPW1hcF9kYmwocHJvYnMsIH4gSFBESSguKVsxXSAlPiUgdW5saXN0KCkpLAogICAgICAgICB1cHBlcj1tYXBfZGJsKHByb2JzLCB+IEhQREkoLilbMl0pICU+JSB1bmxpc3QoKSkKcHJlZHJlc3VsdHMKYGBgCgoKYGBge3J9CnByZWRyZXN1bHRzICU+JSBzZWxlY3QoLWxhbWJkYSwgLXByb2JzKSAlPiUKICBtdXRhdGUodGVtcHM9ZmFjdG9yKHRlbXBzLCBsYWJlbHM9YXMuY2hhcmFjdGVyKHNvcnQodW5pcXVlKGdlcm0uc3RkaSR0ZW1wcykpKSkpICU+JQogIGdncGxvdChhZXMoeD1kYXkseT1tdSxjb2xvcj10ZW1wcyxncm91cD10ZW1wcykpICsKICBnZW9tX2xpbmUoKSAKYGBgCgpBZGQgcmVhbGRhdGE6CgpgYGB7cn0Kc3RkaS5wbG90IDwtIGdlcm0gJT4lIGZpbHRlcihwb3BzPT0iU1RESSIpICU+JSAKICBzZWxlY3QoZGF5LCB0ZW1wcywgY3VtdWxhdGl2ZV9nZXJtLCB0b3RhbF9zZWVkcykgJT4lCiAgbXV0YXRlKHRlbXBzPWFzLmZhY3Rvcih0ZW1wcyksCiAgICAgICAgIHByb3AuZ2VybT1jdW11bGF0aXZlX2dlcm0vdG90YWxfc2VlZHMpCgpwcmVkcmVzdWx0cyAlPiUgc2VsZWN0KC1sYW1iZGEsIC1wcm9icykgJT4lCiAgbXV0YXRlKHRlbXBzPWZhY3Rvcih0ZW1wcywgbGFiZWxzPWFzLmNoYXJhY3Rlcihzb3J0KHVuaXF1ZShnZXJtLnN0ZGkkdGVtcHMpKSkpKSAlPiUKICBnZ3Bsb3QoYWVzKHg9ZGF5LHk9bXUsY29sb3I9dGVtcHMsZ3JvdXA9dGVtcHMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoYWVzKHk9cHJvcC5nZXJtKSwgZGF0YT1zdGRpLnBsb3QpCmBgYAoKUG9vciEKCgo=